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PACS 47 . 27 . E — Turbulent flows, simulation and modeling 
PACS 47 . 35 . Tv - Magnetohydrodynamics in fluids 

Abstract. - The free decay of MHD turbulence at large Reynolds numbers is studied numeri- 
cally using a shell model. We study the statistical properties based on representative sample of 
realisations (128 realisations for each type of initial conditions) over the period of 10 J large-scale 
turnover times. The performed simulations show that the force-free non-helical MHD turbulence 
can demonstrate two different scenarios of evolution in spite of similar initial conditions. Within 
the first scenario, the cross-helicity accumulation is so fast that the energy cascade vanishes before 
significant magnetic energy dissipates. Then the system approaches the state of maximal cross- 
helicity. Within the second scenario, the cascade process continues to remain active until time 
10 4 in units of large-scale turnover time. Then the magnetic field becomes vastly helical due to 
magnetic helicity conservation. Thus the magnetic energy does not dissipate with kinetic energy. 



Introduction. — Free decaying magnetohydrody- 
namic (MHD) turbulence provokes interest for two main 
reasons. First, it raises the possibility of application to the 
physics of the interstellar medium and cosmology in the 
context of evolution of the primordial magnetic field and 
its contribution to the present configuration of the mag- 
netic field in the Universe. Second, the MHD turbulence 
differs from the conventional turbulence of incompressible 
fluids by an extended set of conservation laws, which form 
i the basis of diverse scenarios of free evolution of turbu- 
lent motion. This is a strong fundamental motivation for 
studying this problem. 

Three ideal quadratic invariants are known in 3D in- 
compressible magnetohydrodynamics: the total energy 
E = E v + E b , the cross-helicity H c = (v ■ b) and the 
magnetic helicity H b = {a ■ b), where E" = (|v| 2 /2), 
E b = (|6| 2 /2), v is the velocity filed, b is the magnetic 
field, a is the vector potential (b = V x a). The cross- 
helicity characterises the correlation between the velocity 
and magnetic field pulsations, and magnetic helicity char- 
acterises the correlation between the magnetic field and its 
vector potential. For the non- helical case (H b = H c = 0), 
the decay law for the total energy of the isotropic MHD 
turbulence with equidistributed magnetic and kinetic en- 
ergy (E b (k) w E v (k)) tends to E(t) ~ t' 1 1, which re- 
mains essentially slower than the decay of conventional 
turbulence, for which E(t) ~ t~ 2 . 

The relation between the spectral density of mag- 
netic helicity and the spectral density of magnetic energy 



\H b (k)\ < k~ 1 E b (k) disables the direct spectral transfer 
of magnetic helicity to small scales [5] . Thus the magnetic 
helicity condenses in the largest scale (smallest k) and 
holds there the corresponding part of magnetic energy. 
Then the energy ratio T — E v / E b goes down, and the 
spectral energy flux decreases along with the energy dissi- 
pation rate. Biskamp and Muller p] showed that, for finite 
magnetic helicity, the energy dissipation rate is governed 
by H b and r as e = dE/dt w r x / 2 (l + T)- 3 / 2 E 5 / 2 / H b . By 
making use of numerical result T(t) ~ E(t) they derived 
for r << 1 an asymptotic law E(t) ~ t~ - 5 . Campanelli 
3 specified that the kinetic energy decays as E v {t) ~ t^ 1 
in both nonhelical and helical cases, while the magnetic 
energy follows the kinetic energy in the nonhelical case 
E b {t) ~ E°{t) ~ t' 1 and decays as E b (t) ~ i~ 5 in the 
helical case independently of the initial conditions. Direct 
numerical simulations [1,4 confirmed this phenomenology 
at relatively short time intervals (up to about 10 large- 
scale turnover time), available for computer modelling. 

It is generally believed that for fully developed isotropic 
MHD turbulence there are no reasons for essential correla- 
tion between the pulsations of velocity and magnetic field 
(which means a noticeable level of cross-helicity). How- 
ever, the discovery of highly correlated pulsation of veloc- 
ity and magnetic field in the solar wind [5] awoke inter- 
est in cross-helicity. As far back as 1980 the analysis of 
energy evolution in free-decaying turbulence showed that 
the cross-helicity decays more slowly than energy; there- 
fore the correlation between v and 6, characterised by the 
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normalised cross-helicity C — H c /E, can grow over time 
in free-decaying MHD turbulence [5J. 

The influence of cross-helicity on the forced MHD tur- 
bulence was studied in [7] in the context of the Alfvenic 
turbulence. Using EDQNM it was shown that the system 
tends to a steady state, in which the correlation coeffi- 
cient is much higher than the ratio of helicity to energy 
source injection rates. It has been also determined that the 
energy spectrum of correlated MHD turbulence becomes 
steeper. The role of cross-helicity in stationary forced 
isotropic (not Alfvenic) MHD turbulence has been stud- 
ied by Mizeva et al. [8]. It was shown that the injection of 
cross-helicity suppresses the spectral energy transfer and 
leads to energy accumulation in the turbulent flow. Then 
the spectrum becomes steeper and the intermittency de- 
creases. In the case of stationary forcing with uncontrolled 
injection of cross-helicity, C displays nontrivial behavior 
with long periods of high variability alternating with pe- 
riods of almost constant C [9] . 

Turning back to decaying turbulence, let us note that 
cosmological applications require consideration of very- 
long term evolution of MHD turbulence — the age of 
the Universe calculated in the unit of turnover time for 
the largest galactic turbulent scale (r ~ 10 7 years) gives 
T w 10 4 , which is about 100 times longer than the best 
DNS range. The required time series can be considered 
in the framework of shell models of turbulence. Being an 
ultimate simplification of the original Navier-Stokes equa- 
tions, shell models strictly provide the required conserva- 
tion laws, which are of special interest in the problems 
under discussion. For the first time, the long-time evolu- 
tion of free-decaying MHD turbulence has been considered 
by Antonov et al. [10] . Using the shell model of MHD tur- 
bulence introduced by Frick and Sokoloff [11], they per- 
formed 24 runs of long-time simulation (up to t = 10 
large-scale turnover time) with similar initial conditions. 
The developed Kolmogorov's turbulence with weak mag- 
netic field was taken as an initial condition of the sim- 
ulations. There was no special control of helicities. A 
coherent state with high alignment between the magnetic 
and velocity fields, and essential reduction of the dissipa- 
tion rate, were obtained for most realisations. At the same 
time, there were a few realisations that displayed different 
behaviour, characterised by a low level of cross-helicity. 

In this paper we re-examine the long-time free decay of 
MHD turbulence in the framework of shell models. In con- 
trast to previous work, we use a new (helical) shell model 
of MHD turbulence, introduced in [8], and we study in 
detail the role of magnetic helicity and cross-helicity in 
the evolution scenario. We use a massive computer clus- 
ter which allows us to consider the statistical properties 
based on the representative sample of realisations (128 re- 
alisations for each type of initial conditions) . 

Shell model for helical MHD turbulence. Shell 
models describe the dynamics of fully-developed MHD tur- 
bulence through a set of complex variables U n , B n , which 



characterise the amplitudes of velocity and magnetic field 
pulsations in the shell of wave number k n < \k\ < k n+Xl 
where k n = A" ( A is the shell width in a logarithmic scale) . 
The total energy in terms of shell model is obviously the 
sum of kinetic and magnetic energy of individual shells 
E = J2(\U n \ 2 + \B n \ 2 )/2. The cross-helicity is defined in a 
similar way as H c = £(E/„-B* + B n U*)/2. The definition 
of magnetic helicity is not so evident. The most popular 
MHD shell models [TT ] fT2 ] introduce the magnetic helicity 
as H b — J2(~l) n \B n \ 2 / k n , which associates the magnetic 
energy of a given shell with a positive or negative mag- 
netic helicity. Then the non-helical state can be obtained 
only through the balance of magnetic energy in the even 
and odd shells. Following the idea used for the hydrody- 
namical helicity in [13j , we define the magnetic helicity as 
H b = E fc n 1 ((- B n) 2 - B n)/ 2 - Thi s definition allows us to 
get helicity of any sign in any shell. So we use the model 
introduced in [5]: 

ckU n =ik n (K(U,U)-K(B,B))-^, (1) 

Re 

d t B n = ik n (A n (U, B) - A n (B, Uj) - (2) 

Km 

where the nonlinear terms are written as 
A n (X, Y) — \ 2 {X n+x Y n+x + X* +1 Y* +1 ) — X r n _ x Y n 

+X n - 1 Y n - 1 + X*_ X Y*_ X - \ 2 (X'n +1 Y n + X n Y^ +x ) 
+i\(2X*Y£ +1 + X r n _ x Y^_ x — X^_ X Y^_ X ), 

a star means complex conjugation, and superscripts r, i are 
real and imaginary parts. Re and Rm are the kinetic and 
magnetic Reynolds numbers. In the limit Re,Rm — > oo 
these eqs. {J) and ([2]) conserve the total energy, the cross- 
helicity, and the magnetic helicity. Time is measured in 
dimensionless units equal to the eddy turnover time on 
the shell n = (fco = 1), which corresponds to the largest 
scale of the system. For the case of hydrodynamics, L'vov 
et al. [T3] suggested the parameter A equal to the golden 
number (1 + V5)/2 for optimal spectral resolution which 
is also successfully applicable for MHD turbulence [15] . 

Numerical results. — Equations ([1]) and were 
integrated for < n < 40 up to the time t = 10 5 , when 
any turbulent transfer has been finished and only pure ex- 
ponential energy decay remains in the largest scales. In 
all simulations, A = 1.618 and Re = Rm = 10 5 . First we 
consider the decay of fully-developed turbulence with van- 
ishing helicities H c and H h . Initial values of shell variables 
are Uq = — ^^/2 + 8, B a = \J~2 + <5, where S is a random 
complex additive with real and imaginary parts in a range 
[-10~ 4 : 10- 4 ], and U„ = B n = for other shells (n > 0). 
It corresponds to the initial state with E u « 1, E b ss 1 and 
a small quantity of helicities |if c | < 10 4 , | < 10~ 4 . 

The inertial range with a Kolmogorov "-5/3" power law 
is formed at about the turnover time of the vortex of maxi- 
mal scale (t w 1). In fig. [T] we present the energy spectrum 
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Fig. 1: Evolution of the compensated power spectrum, aver- 
aged on 128 realisations: t = 1, 10, 100, 10 3 , 10 4 (from top to 
bottom). 



averaged on all 128 realisations at different stages of evo- 
lution. The spectrum is compensated by fc 2 / 3 which leads 
to a horizontal line for the Kolmogorov spectral law. The 
inertial range practically disappears at the time t w 10 3 , 
although the sagging in the middle part of the spectrum 
indicates that there are some realisations in which the en- 
ergy is present at relatively small scales. 

Figure [U[a) shows the evolution of the total energy in 
the entire totality of realisations. We see that similar ini- 
tial conditions lead to the different scenarios in the evo- 
lution of the system. At the early stage of development 
(on the periods from several units to several ten) the bun- 
dle of trajectories remains sufficiently dense and is limited 
by the power laws E(t) ~ t^ 1 and E(t) ~ t~ x / 2 , shown 
in the figure by dashed lines. Note that these are two 
power laws suggested for the decay of non-helical and he- 
lical MHD turbulence (in sense of magnetic helicity) [3]. 
At the next stage (50 < t < 1000), the bundle of trajecto- 
ries is bounded from below by a power law as before (but 
more steep, like E(t) ~ t~ 2 ), while from above separate 
trajectories leave the bundle practically horizontal, which 
indicates the vanishing cascade of energy and transition 
to the exponential dissipation of the energy. 

The variety of scenarios of the evolution confirms 
fig- Efb) j m which we show the evolution of the normalised 
cross-helicity C. The quantity C = H c / E characterizes 
the part of energy, concentrated in the correlated pulsa- 
tions of velocity and magnetic field. The limit C — > ±1 
corresponds to a completely correlated state (U n = B n ), in 
which the nonlinear energy transfer is blocked. Figure^b) 
shows that the main part of the trajectories reaches this 
state in the range of time 100 < t < 1000, but there are 
some trajectories (about 10-15%) for which the evolution 
of C stops at some arbitrary level. 

Figure OHa) shows that the cross-helicity H c can be gen- 
erated, in contrast to energy, which can decay only. The 
active cross-helicity production is mostly observed until 
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Fig. 2: Evolution of total energy (a) and normalised cross- 
helicity (b). 128 realisations with similar initial conditions. 



t w 10 2 . In some realisations the cross-helicity reaches a 
level \H C \ k, 0.2. The source of the cross-helicity in the 
force-free evolution can be the dissipation term only. In 
fig- Etb) we present the evolution of the spectrum of cross- 
helicity for one realisation (namely, we took the realisation 
that corresponds to the lowest trajectory in fig-EJa)). One 
can see that most intensive production of the H c happened 
at the dissipation scale. Then nonlinear terms transfer H c 
to the largest scale like an inverse cascade. Note that in 
contrast to the curves in fig.[IJ the lower curves in fig.[3][b) 
correspond to earlier moments of time. 

The cross-helicity, which can be produced in the small- 
est scales by dissipation, is transported towards large 
scales through the inertial range. Then the spectral energy 
flux is considerably reduced depending on the level of ac- 
cumulated cross-helicity in a given realisation. This leads 
to the blocking of the turbulent energy cascade, which 
occurred in different realisations at substantially different 
moments in time and with substantially different values 
of the remaining total energy (see fig. [2^). However, the 
growth of the normalised cross-helicity does not result in 
the completely correlated state C = ±1 for all realisa- 
tions (see fig. EJd) . Some realisations continue to stay at 
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Fig. 3: Evolution of the cross-helicity H c for all realisations (a) 
and spectral distribution for a single realisation at time t=0.8, 
0.9, 1, 5 (lines from bottom to top) (b). 



values — 1 < C < 1. These realisations develop a highly 
helical magnetic field. The magnetic helicity (which is 
weak at the initial state) does not cascade to small scales 
and practically does not dissipate. This has been demon- 
strated through various studies starting from [2]. If the 
energy transfer (and dissipation) is not blocked by the 
cross-helicity, only the magnetic field with maximal helic- 
ity survives at the late stage of the evolution. 

The tendencies of correlation evolution are illustrated 
fig. [4j in which the distribution of realisations is shown 
on the (C 6 ,C) plane at different moments in time (here 
C b = H b /(k E b ) ). All realisations start from the origin 
of the graph (at t = 0). At the first stage (up to t « 100), 
points scatter along the vertical line, showing the rapid 
increase of |C|. Several realisations (« 10%) deviate from 
the axis C b — 0. At the time t = 10 3 , this set of points 
forms a cloud in the center of the plane, while the rest are 
concentrated on the lines C = ±1. At the late stage (t = 
10 4 ), practically all points are on the lines C — ±1 or C b — 
±1. The first case (C — ±1) means that magnetic and 
velocity fields are completely correlated, while the second 
case (C b = ±1) means that only the helical magnetic field 
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Fig. 4: normalised cross-helicity C vs normalised magnetic he- 
licity C for different realisations at t = 10 2 (crosses), t = 10 3 
(open circles) and t = 10 4 (full circles). 



remains at the largest scale. 

Figure [5] shows the relation between the kinetic and 
magnetic energies for all runs at different times. One can 
see that both energies scatter randomly from 0.01 to 0.1 at 
early stage (t — 100). At time t — 10 3 , kinetic and mag- 
netic fields approach the equipartition state (U n = B n ) 
and the points mainly lie close to the line E b = E u . At 
late stage (t = 10 ) some points move away from this 
line. These points correspond to realisations that follow 
the second scenario, in which the kinetic energy continues 
to decay while the magnetic one, being helical, does not 
change (C fc = ±1). 

Discussion and Conclusions. The performed 
simulations show that the force-free non-helical (H b ss 
0, H c 0) MHD turbulence can demonstrate fundamen- 
tally different ways of evolution in spite of similar ini- 
tial conditions. We distinguish two scenarios of evolution. 
Within the first scenario, the cross-helicity accumulation is 
so fast that the energy cascade vanishes before significant 
magnetic energy dissipates. Then the system comes to a 
state with C = ±1 and the value of C b depends on the rest 
of the magnetic energy. Within the second scenario, the 
cascade process remains active until the late time t ~ 10 4 , 
when the magnetic field becomes vastly helical and later 
magnetic energy does not dissipate with kinetic energy. 
Then the system comes to the final state with C b — ±1 
and an arbitrary value of C. 

Comparison of the results with a direct numerical simu- 
lations is hardly possible, because in high resolution DNS 
the reachable time is typically less than 100 dimensionless 
units. Shell model simulations have shown that the dif- 
ference in evolution scenario is well pronounced after 10 3 
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Fig. 5: Magnetic energy i? b vs kinetic energy i? u for different 
realisations at t = fO 2 (crosses), f = fO 3 (open circles) and 
t = fO 4 (full circles). 

units of time (see fig. E]) . 

The probability of appearance of both scenarios depends 
on the presence of helicities H c and H b in the initial distri- 
bution. We performed simulations with initial conditions 
E u (k ) — E b (k ) — 1 and some given value of C or C b 
in all 128 realisations. It was found that the presence of 
even a small portion of magnetic helicity substantially in- 
creases the probability of the second scenario: C b = 0.01 
at t — leads the second scenario being developed in half 
of the cases, and initial value C b = 0.05 leads to the likeli- 
hood that practically all realisations will follow the second 
scenario. The presence of cross-helicity at the beginning 
of evolution is not so crucial (because the turbulence it- 
self generates the cross-helicity). C = 0.01 at t — does 
not change the distribution of trajectories at all. Only 
C ~ 0.1 is enough to exclude the appearance of the sec- 
ond scenario. 

We stress the point that the dynamics of helicities has 
a signifiable influence on the evolution of free-decaying 
MHD turbulence. Our understanding is based on the idea 
of the inverse cascade of magnetic helicity [2] which has 
been confirmed by numerous simulation [4| l 16< - fl9] . Alex- 
akis et al |19] reported that smaller-in-amplitude direct 
cascade is observed from the largest scale to small scales. 
This may be a result of insufficient inertial range resolu- 
tion, in which forcing and dissipation scales are not well 
separated. In shell models, the Reynolds number is always 
quite large (10 5 or more) and the direct cascade of mag- 
netic helicity was never observed. From the other perspec- 
tive, a shortcoming of shell models is ignoring the nonlocal 
(in scales) interactions, due to which weak direct cascade 
of magnetic helicity may occur. This shortcoming can be 
overcome in the framework of shell models - a detailed 
study of the shell-to-shell interactions [3U] gives a base on 
which to build a nonlocal MHD shell model [21] . 

Finally, we note that the relationship \H b (k)\ < 
k~ 1 E b (k) indicates that the character of evolution can be 
changed substantially if the magnetic helicity can move to 



U 0.0 




1 10 100 1000 if) 4 10 5 



Fig. 6: Evolution of normalised cross-helicity in a sample of 128 
realisations in the system with extended large-scale spectral 
range (k min =0.1). 

scales larger than the scale of maximal energy at t = 0. 
To examine this case, we extended the range of scales into 
the red part of the spectrum, simulating eqs. (JlJ and © 
for — 5 < n < 35. This means that the maximal scale ac- 
cessible for the turbulence becomes about 10 times larger 
than the scale in which the energy is set at the initial state 
(kmin ~ O.lfco). The result of the simulation is shown in 
fig. El In contrast to fig. EJb), there are no more real- 
isations that avoid the completely correlated final state 
(C = ±1). The inverse cascade of magnetic helicity leads 
to the fact that an increasingly smaller part of magnetic 
energy might be blocked at the largest scales, being ex- 
cluded from the direct energy cascade to small scales. 
Then only the correlated part of velocity and magnetic 
fields survives at the late state of evolution. 
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